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Abstract 

We study numerically and theoretically (on a heuristic level) the time evolution 
of a gas confined to a cube of size L 3 divided into two parts by a piston with 
mass Ml ~ L 2 which can only move in the x-direction. Starting with a uniform 
"double-peaked" (non Maxwellian) distribution of the gas and a stationary piston, 
we find that (a) after an initial quiescent period the system becomes unstable and 
the piston performs a damped oscillatory motion, and (b) there is a thermalization 
of the system leading to a Maxwellian distribution of the gas velocities. The time 
of the onset of the instability appears to grow like L log L while the relaxation time 
to the Maxwellian grows like L 7 /2. 

1 Introduction 

The time evolution of a gas filled container divided by a massive piston is an old problem 
which has attracted much attention recently from both a conceptual and computational 
point of view, [jO], O |KBM|| . While we do not believe that there is any "paradox" 



associated with this problem, the basic mechanism of energy transport across the piston 
from the hot side to the cold side which are at equal pressure but different temperatures, 
was already described for an idealized version, by one of us in 1959 \fA \, many intriguing 



questions remain. Not surprisingly the microscopic motion of the piston, and thus the 
time evolution of the system from an initial, far from equilibrium state, to the final 
true equilibrium state is not easy to compute analytically or even numerically for large 
systems [[KB1\1|| . It is therefore interesting to consider the greatly simplified case when 



the gas particles only interact via collisions with the massive piston [|LT| , |GF| , |GP|| . This 



department of Mathematics, University of Alabama at Birmingham, Alabama 35294 
2 Department of Mathematics, Rutgers University, New Jersey 08854 
3 Current address: Institute for Advanced Study, Princeton, NJ 08540 



1 



serves, among other things, as a model for the approach to equilibrium in such a "weakly" 
interacting macroscopic system. 

Here we carry out studies on this system under the "simplest" initial conditions. 
The positions and velocities of gas particles of unit mass, in an insulated box of size 
L 3 , are picked, at t — 0, from a Poisson process with some density function, i.e. they 
are assumed to be independent random variables with a given double peaked velocity 
distribution p(v) = p(—v) (see (|2.4j )) and a constant spatial density throughout the box. 
The piston, which acts as a single particle with area L 2 and mass Mi ~ L 2 , is then 
released at the position X(0) = L/2 with zero velocity V^(0) = 0. 

It might be thought that the random fluctuations of the piston's velocity due to the 
initial randomness in the positions and velocities of the gas particles will be relatively 
small when L is large (the number of gas particles growing like L 3 ) and vanish in an 
appropriate hydrodynamic scaling limit (as L — > oo) of time, space and piston mass. 
In fact, some rigorous results in this direction were obtained in [ |LP!j| ] and |CL!a| ] for a 



class of initial conditions more general than those considered here. It is proven in those 
papers that for those initial distributions the dynamics of the piston, in the hydrodynamic 
scaling limit, is governed by deterministic equations for as long as no particle collides 
with the piston more than twice. After that period, however, we could not control the 
random fluctuations in the particle configuration anymore. It is clear, though, that for 
an initial state which is spatially uniform and has the same velocity distribution on 
both sides of the piston, the deterministic macroscopic evolution would predict that the 
position of the piston remains stationary forever. On the other hand, when the initial 
velocity distribution of the particles is not Maxwellian, the system is not at thermal 
equilibrium, and we expect it to somehow evolve toward equilibrium for almost any 
initial configuration (with respect to the Liouville measure on the energy surface). This 
is indeed what we found for every initial distribution. The path to the Maxwellin taken 
by the system turned out however to be quite sensitive to the initial distribution. 

Here we present the results of detailed numerical simulations for one particular initial 
density p(v) that vanishes unless 0.5 < \v\ < 1. We found that the quiescent state 
becomes unstable after a few (5-10) recollisions of each gas particle with the piston. 
The system then develops large (on a macroscopic scale) oscillatory motions which last 
for a long time. They look like smooth harmonic oscillations with an (almost) constant 
amplitude and piston speed comparable to that of the gas particles. Eventually, though, 
the oscillations dampen and vanish, and the system approaches a stable equilibrium state 
with constant gas density and Maxwellian velocity distribution0. 

Our main conclusions from the numerical investigations, which we also "derive" 



1 The latter is of course expected on general grounds for, as pointed out by Boltzmann: "[the Maxwell 
distribution] is characterized by the fact that by far the largest number of possible velocity distributions 
have the characteristic properties of the Maxwell distribution, and compared to these there are only 
a relatively small number of possible distributions that deviate significantly form Maxwell's. Whereas 
Zermelo says that the number of states that finally lead to the Maxwellian state is small compared to all 
possible states, I assert on the contrary that by far the largest number of possible states are "Maxwellian" 



and that the number that deviate from the Maxwellian state is vanishingly small", B, see also [L2] 
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heuristically, is that 1) while the dynamics are non chaotic in the technical sense of there 
being no positive Lyapunov exponents, there are for typical initial conditions! enough 
interactions to bring the system to equilibrium along "interesting" nontrivial pathways, 
and 2) that the time of onset of the instability grows with L as LlogL. This means that 
even on the hydrodynamical scale r = t/L, the deterministic behavior in which nothing 
happens on the spatial scale, y = x/L, would remain valid when L — ► oo. In terms 
of r however this onset time grows only like logL so that in "practice" this behavior 
extends to macroscopic systems. We believe that this instability is related to the fact 
that the deterministic solutions are unstable for the kind of p(v) we consider here. For 
other initial p(v), for which the deterministic solution is stable, we do not get such large 
oscillations but, as expected, we still get an approach to Maxwellian, see Section 7. 



2 Description of the model 

Consider a cubical container = [0, L] x [0, L] x [0, L] filled with an ideal gas consisting 
of N particles. The container is divided into two parts by a wall (piston) orthogonal to 
the x axis. At time t = the wall is released and then it can move freely without friction 
inside A^ along the x-axis, under the action of elastic collisions with the gas particles, 
each of which has the same fixed mass m. Since the piston's area is L 2 , we assume its 
mass Ml to be proportional to L 2 and given by Ml = hmL 2 with a fixed constant b > 
(we set m — 1 and h = 2 in our numerical simulations). 

Since the components of the particle velocities perpendicular to the x-axis play no role 
in the dynamics of the piston, we may assume that each particle has only one component 
of velocity, v, directed along the x-axis. Hence, each gas particle can be specified by 
a pair (x«, t>j), where i = 1, . . . , N. We shall take the total number of particles A" and 
the total kinetic energy of the system proportional to L 3 ; and we are interested in the 
behavior of the system for L 3> 1. 

The initial configuration of gas particles and their velocities {{xi,Vi)}f =1 is selected 
at random with statistics given by a (two-dimensional) Poisson process on the x, v plane 
with a given density, p(x,v) = L 2 p(v). More precisely, for any domain D C [0,L] x H 1 
the number of particles in D, i.e. N D = : (x^Vi) G D}, has a Poisson distribution 
with parameter 

Xd = L 2 / p(v)dxdv 
Jd 

That is, for each k > 0, 

\ k 

P(N D = k) = -fe- x ° 
k\ 

2 Typical here means of probability converging to one (presumably, exponentially fast) as N — > oo 
with respect to the invariant Liouville measure. There are clearly some exceptional initial configurations 
(e.g., those with a complete symmetry about x = L/2 with opposite velocities) when V(t) = at all 
times and so the speed of gas particles never changes. We expect (for the reasons given by Boltzman 
Q) all such initial states to be unstable with respect to generic small perturbations. 



3 



The numbers of particles in nonoverlapping domains are independent. The function 
p(v) = p(—v) takes values of order one, and the factor of L 2 is simply the cross-sectional 
area of the container Ax,. The piston is initially at rest at the midpoint, X(0) — L/2, 
V(0) = and the position and velocity at time t are denoted by X, < X < L, and V. 

We note that the total number of gas particles A" in the container is random. So are 
the numbers of particles to the left and to the right of the piston, call them N_ and A^ + , 
respectively (N = N_ + N + ). The initial total kinetic energy of the system, \ J2 mv h * s 
also a random variable. Once chosen, the values of AL , N + and the total kinetic energy 
E of gas plus piston are (presumably, the only) integrals of motion. 

The gas particles and the piston move freely (with constant velocity) between elastic 
collisions of particles with the piston and the walls. When a particle collides with a wall 
at x = or x = L, its velocity simply reverses. If a particle with velocity v hits the 
piston whose velocity is V, then their velocities after the collision, call them v' and V, 
respectively, are given by 

V'=(l-e)V + ev (2.1) 
v' = -(1 - e)v + (2 - e)V (2.2) 

where 

2m ( 2 \ . . 

(2.3) 



M + m \bL 2 + l, 

In order to avoid recollisions of gas particles with the piston, for at least some initial 
period of time, we impose a velocity cutoff! 

p(v) = if \v\ > w max or \v\ < v min (2.4) 



with some < t> m i n < f max < oo, cf. [ |ULS| , |LPS|| . Under these conditions, there will 



be an initial time interval of length O(L) during which each particle colliding with the 
piston will have to travel to the wall, bounce off it, and then travel back to the piston 
before it can hit it again. Therefore, during that interval, there will be no recollisions 
of any gas particle with the piston. We call such an interval of time the zero-recollision 
interval. Likewise, if the piston remains slow enough after recollisions occur, there will 
be another interval of time of order L during which each gas particle experiences at most 
one recollision with the piston. We call it the one-recollision interval. 

The hydrodynamic limit is now obtained by rescaling space, y = x/L, < y < 1, and 
time, t = t/L. In the new space-time coordinates y, r the zero-recollision interval and 
the one-recollision interval will be of order one, and we will denote them by (0, To) and 
(r ,ri) respectively. As already mentioned, we prove in ||CLS|1 , for a more general class 



of initial densities p(x,v) that during the interval (0, n) the functions Yl(t) and Wl(t) 
Y L (r) = X(tL)/L, W l (t) = dY L /dr = V(rL) 



3 As explained in Section 7, the choice of p(v) plays an important role in the later time evolution of 
the system. 
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converge uniformly in probability, as L — > oo, to some deterministic functions Y(r) and 
W(t); W satisfies a certain algebraic equation and dY/dr = W. In the case considered 
here, when the density function p(x, v) is independent of x and symmetric in v, i.e. 
p(v) = p(—v) for all x, v, then the hydrodynamic evolution is trivial: Y(t) = 0.5 and 
W(r) = for all r > 0. Hence, Yl{t) — > 0.5 and Wl{t) — > as L — > oo, uniformly for 
all r e (0,ri). 

In this paper we investigate numerically what happens to Y L {r) and Wl(t) beyond the 
interval (0, 7~i). More precisely, for how long does the stochastic trajectory {Y L (r), Wl(t)} 
remain close to the deterministic one {W(t),Y(t)}7 And what does it look like in the 
long run? 



3 Numerical results 

In the computer simulations described here and in Sections 4-6, we set 

*(«) = { I elsethei^' " 1 (31) 

so f m i n = 0.5 and t> max = 1. The x and v coordinates of all the particles are then 
independent random variables uniformly distributed in their ranges < x < L and 
Vmm < \v\ < fmax- Our computer program first selects N according to the Poisson law 
with mean L 3 , and then generates all (xi,Vi), 1 < i < N, independently according to 



their uniform distributions. We used the random number generator described in |[MN|| . 
The parameter L changed in our simulations from L = 30 to L = 300. For L = 300 the 
system contains ~L 3 = 27, 000, 000 particles. 

Once the initial data is generated randomly, the program computes the dynamics by 
using the elastic collision rules ( |2.1|) , (|2.2| ). All calculations were performed in double 
precision, with coordinates and velocities of all particles stored and computed individ- 
ually. We note that memory requirements alone can be enormous - the program needs 
over 430Mb RAM in order to run the model with L = 300. Also, the task of determining 
which particle is to collide with the piston next is nontrivial. To avoid a long search 
through the entire set of iV particles after each collision, we assembled a smaller group of 
particles located in a vicinity of the piston, where the search is performed, and updated 
this group periodically, as time goes on. The calculations were done in C++ and HPF 
(High Performance Fortran) on a Dell Power Edge machine with Dual 733MHz proces- 
sors at the University of Alabama in Birmingham. Our code is available on the web page 
referred to in the next section. 

Figure 1 presents a typical trajectory of the piston. Here L = 100. The position and 
time are measured in hydrodynamic variables Y = X/L, < Y < 1, and r = t/L. 

Initially, the piston barely moves about its equilibrium position: recall that the hy- 
drodynamic trajectory of the piston is Y(t) = 0.5 for all r > and that this holds 



exactly for r < 2 as L — > oo, ||CJLS |. Then, at times r between 3 and 5, the random 
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vibrations of the piston grow and become quite visible on the y-scale, but for a short 
while they look random. After that the piston starts travelling back and forth along the 
y axis, making excursions farther and farther away from the equilibrium point y = 0.5. 
Very soon, at r = r max « 8, the swinging motion of the piston reaches its maximum, 
(Ay) max = max \Y(t) — 0.5 1 ~ 0.1. Then the oscillations of the piston dampen in size 
and seem to stabilize at an amplitude A « 0.04. At the same time the trajectory of the 
piston smoothes out and enters an oscillatory mode with a period r per « 1.63. 

The velocity of the piston W(t) follows a similar pattern. Its random fluctuations 
grow after r ~ 2.5, then it starts swinging up and down, reaches the maximum value of 
W max = max |W(r)| ~ 0.4 at time r m 9.5. After that the oscillations of W(r) dampen 
and seem to stabilize. 

Note that the graph of the function Y(r) looks much smoother than that of W(t), as 
would be expected from the fact that Y(r) is the integral of W(r). Interestingly, both 
functions Y and W smooth out as time goes on. 

This cycle of the gas motion in the container between the walls and the piston con- 
tinues for a long time with the same period r pcr ~ 1.63, independent of L, but the 
amplitudes of both Y(r) and W(r) are slowly decreasing, see Fig. 3. 

The oscillations of the piston with decaying amplitude can be described, in the interval 
30 < r < 1000, approximately by 

Ki(r) ~ Ae- x{T - 20) smu(r - a) (3.2) 

with A = 0.046 and some constant A > 0. Correspondingly, W\{t) = dYi/dr in the same 
interval 30 < r < 1000 is 

W x {t) ~ -\Y 1 + Ae~ KT - 20) ucosuj{r-a) 

= Ae~ x(T ~ 20) [-\smuj(T - a) + u cosuj(t - a)} 

= A ie - x ^ 20 hmuj{T - (3) (3.3) 

with A\ = A-JuS 1 + A 2 and some (3 related to a. 

To check how well our formula (|3.2|) agrees with the data, we computed the amplitude 
A(t) as a function of time r, by fitting a sine function Y q (t) = Asib.uj(t — a) "locally", 
on the interval (r — 5, r + 5) for each r. Fig. 4 shows A{r) on the logarithmic scale, 
which looks almost linear on the interval 30 < r < 800. (After that, Y{r) becomes quite 
unstable, with already small amplitude A(r) decreasing abruptly, possibly due to the 
interference from random fluctuations, so we left that part out.) 

We used the least squares fit to estimate A = 0.00264 for the run shown on Figs. 1-4. 
Since A is small, the oscillations indeed die out very slowly. The "half-life" time (the time 
it takes to reduce the amplitude by a factor of two) is tu^ = A -1 In 2 w 263. Note that 
over the time interval 30 < r < 800 we only observed the reduction of the amplitude by 
a factor of about 10, and the periodic oscillations were still visible on the plot at r > 900. 
Also, A and hence T1/2 depend on the system size L, see below. 
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4 Dependence of the Time Evolution on System 
Size 



Many of the characteristics of the piston trajectory described above ((AY) max , W ma x, A 
Tper) appear to be independent of L. The following tabled presents computed values of 
all these parameters for different L's. 



L 


(^^)max 


w max 


A 


Tper 


30 


0.114 


0.40 


0.042 


1.70 


50 


0.121 


0.39 


0.045 


1.65 


80 


0.105 


0.37 


0.042 


1.65 


100 


0.100 


0.34 


0.041 


1.64 


120 


0.122 


0.39 


0.041 


1.62 


150 


0.111 


0.37 


0.045 


1.62 


200 


0.102 


0.37 


0.037 


1.62 


250 


0.100 


0.41 


0.042 


1.64 


300 


0.122 


0.42 


0.045 


1.65 



Table 1. Principal characteristics of the piston dynamics. 

In each case, we averaged over several experimentally generated trajectories (for L = 
250 and 300, just one trajectory was used). 

There are however other quantities such as r max , n/2, and the related A, which depend 
in a systematic way on L. In particular, we estimated numerically that r ly / 2 ~ L , hence 
A ~ L~ L3 , see Fig. 5. We will denote A = A^ and discuss it further in Section 6. We also 
noticed that in some runs with larger L's (such as L = 150 and L = 200) the exponent 
A changes with time, it is higher when r < 100 and lower when r > 100. 

But most importantly, the time of the largest oscillations r max and the related time 
of the onset of the instability r c , see below, seem to slowly grow with L, very likely as 
logL. To understand this fact, we looked into the mechanism of the build-up of random 
fluctuations of the piston position and velocity, which eventually result in their large 
nearly harmonic oscillations. To this end we plotted the histogram of the (empirical) 
density of gas particles in the y, v plane at various times < r < 30, see samples in 
Fig. 6. The initial density (at time zero) is almost uniform over the domain < x < L 
and f m i n < \v\ < f m ax (variations in the initial configuration always exist, because it is 
generated randomly). Then, for < r < 1, the piston experiences random collisions with 

4 The numerical data here are given primarily for demonstrating the typical scale of oscillations. They 
are not meant to be estimates of physical parameters, so we do not provide error bars. For L < 200, 
where more than one experimental trajectory was generated, we have estimated that the accuracy of 
our numerical values in Table 1 is at least within 10%. 
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particles and acquires a speed of order M L = 0(1/L), see [|L~T| , po| , [DGL|| . These small 



fluctuations of the piston velocity result in changes of the velocities of the particles which 
leave the piston after collision. Thus the outgoing particles on the right hand side of the 
piston have velocities in the interval (v m i n + 2W(r), v max + 2W(r)) while those on the left 
hand side of the piston have velocities in the interval (— f m in + 2W(r), — v max + 2W(r)). 
Hence, the region in the y, v plane where the density of the particles is positive is no 
longer a rectangle with straight sides, now its boundaries are curves whose shape nearly 
repeats the graph of the randomly evolving piston velocity W(r). While the variations 
of 0(1/ L) of these boundary curves may seem small, it is crucial that on opposite sides 
of the piston they go in opposite directions. Indeed, when W(r) > 0, then the outgoing 
particles on the right hand side accelerate and those on the left hand side slow down. 
When W(t) < the opposite happens. 

Next, the particles that have collided with the piston travel to the wall and come 
back to the piston. Now their densities are less regular than they were initially - the 
regions in the x, v plane where the density is positive, are curvilinear domains. When 
they hit the piston, they shake it back and forth more forcefully than before, because the 
velocities of the incoming particles on the opposite sides of the piston are now negatively 
correlated. When particles on the right hand side are fast, those on the left hand side 
are slow, and vice versa. This produces a resonance-type effect destabilizing the piston 
dramatically and the velocity of the piston W(t) experiences larger fluctuations than 
before. The velocities of the newly outgoing particles will again go up and down in 
opposite direction, on a greater scale than before. 

As time goes on, the above phenomenon repeats over and over, with larger and 
larger fluctuations of the gas and piston velocities, until the distribution of gas particles 
completely breaks down. At times r ~ 10, two large clusters of particles are formed, 
one on each side of the piston. When one cluster bombards the piston, the other moves 
away from it and hits the wall, then they exchange their roles. The clusters have sizes 
of about 0.3-0.5 in the y direction and the particle velocities range from about 0.2 to 
just over 1. The average velocity is about 0.5-0.6 and so the clusters hammer the piston 
periodically with period 1.6-2.0, which is close to the experimentally determined period 
of piston oscillations, see above. 

Fig. 6 shows six snapshots of the empirical density of gas particles taken at different 
times. At r = the gas fills (almost uniformly) two rectangles {(y, v) : 0.5 < \v\ < 1, < 
y < 1}. At t = 2.3 one can see some ripples on the boundaries of these rectangles. At 
time r = 4.2 the irregularities grow and at r = 5.9 the rectangular shape is broken down. 
Two large clusters of particles are formed, both appear in the upper half-plane v > 0, i.e. 
at that time both clusters move to the right (one toward the piston, the other away from 
it). Later the density undergoes strange formations (r = 7.4) but eventually smoothes 
out and enters a slow process of convergence to Maxwellian (r = 18.6) described below. 
Note a narrow white line around v = 0, meaning the total lack of very slow particles at 
time r = 18.6. 
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A longer sequence of snapshots at times < r < 30 is posted on the web page 



www.math.uab.edu/chernov/piston/pictures/piston.html 

It gives a spectacular view of the entire system evolution. 

The above analysis may suggest that the fluctuations of the piston velocity roughly 
increase by a constant factor during each time interval of length one. Indeed, initial 
random fluctuations W a ~ 0(1/ L) result in additional changes of velocities of outgoing 
particles by 2W a . When those particles come back to the piston (in time At « 1), 
they kick its velocity to the level of 2W a . Then the newly outgoing particles acquire 
an additional velocity 4W a , etc. Over each time interval of length one the fluctuations 
double in size. This is an obvious oversimplification of the real dynamics, but it leads to 
a reasonable conjecture 

W a (r) » ^ (4.4) 

where W a (r) are typical fluctuations of the piston velocity at time r and C, R > are 
constants. 

We tested the above formula numerically as follows. Let W c > be some preset 
critical value of the piston speed and r c = inf{r > : |VK(r)| > W c } the (random) time 
when W c is first reached. This time plays the role of the "onset" of large fluctuations of 
the piston velocity. One would expect, based on (|4.4|) that 



T c &\n(W c L/C)/lnR (4.5) 

i.e. r c grows as InL when L increases. 

We found r c experimentally for W c = 0.1 and W c = 0.15 and checked that (|4.5| ) 
agreed well with the data, see Fig. 7. By the least squares fit we estimated C = 0.45 and 
R = 1.6. 



5 Approach to equilibrium 

We examined the convergence of the velocity distribution of gas particles to a Maxwellian. 
At any given time r > 0, let 

F r (u) = #{* : Vi < u}/N 

be the empirical (cumulative) distribution function of particle velocities. At equilibrium, 
it should be close to the normal distribution function $(x). As a measure of their 
closeness, we used the supremum of the difference 

D T = sup \F T (u) - $(it)| 

— oo<«<oo 

Initially, Do ~ 0.245 for our choice of p(v). One can expect that D T converges to as r 
grows when iV is large. In fact, if the velocities Vi were independent and had Maxwellian 
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distribution (which it would be in true statistical equilibrium), then D T would be of 
order 0(l/y/N), and the product D T \^N would have a certain limit distribution, see the 
theory of the Kolmogorov-Smirnov statistical test ||hu|| . In particular, it is known that 



for a Maxwellian the probability P(D T y/N > 1) « 0.2. Based on this, we define the time 
of convergence to equilibrium by 



r, 



eq 



inf{r > : D T VN < 1} (5.1) 



Here the constant 1 as a critical value is chosen arbitrarily. We estimated r eq for various 
L's and found that r eq « aL b with some constants a,b > 0. By a least squares fit we 
found a = 0.18 and b = 2.47, see Fig. 8. (Note that the accuracy of our experimental 
data seems to increase with L, as the points are getting closer to each other and to the 
fitting line for larger L's on Fig. 8.) 

The plot of the product S = D T yN versus r is given on Fig. 9 (for a particular run 
with L = 40). It shows that, after an initial sharp drop for < r < 20, the statistic S 
decreases exponentially in r. Another commonly used statistic to measure closeness to 
Maxwellian is 

S' = 3- 



M 4 
Ml 



where Mi and M4 are the second and the fourth sample moments of the empirical veloc- 
ity distribution, respectively. Fig. 9 shows that S' converges to zero in a similar manner 
(for the same run with L = 40). 

Theoretical Considerations 

As we noted in Section 0, the total energy E and the numbers of particles in the left 
and right compartments (iV_ and iV + ) are integrals of motion. With these quantities 
fixed, the model can be reduced to a billiard system in a high-dimensional polyhedron 
by standard techniques, as we show next. 

Let {xi}, i = 1, . . . , N + , denote the x-coordinates of the particles to the right of the 
piston, and {xi}, i = — 1, . . . , — AL, those to the left of it (ordered arbitrarily). Put 
Xq = XyfM, where X is the coordinate of the piston and M is its mass. Then the 
configuration space of the system (in the coordinates Xj, — iV_ < i < N + ) is a polyhedron 
Q C !R Ar+1 (recall that N = N_ + N + ) defined by inequalities 



< X-N-, ■ ■ ■ ,X-i < xq/V M < x\, . . . ,xn + < L 

It is known that the dynamics of our mechanical system ( "gas+piston" ) corresponds to 
the billiard dynamics in Q, see flCFSfl . That is, the configuration point q e Q moves 



freely and experiences specular reflections at the boundary dQ. The velocity vector 

p = q = {v- N _,. . • VVM, Vi, . . .,v N+ } 

has constant length, since ||p|| 2 = 2E = const. Therefore, the phase space of the billiard 

'p where 



system is M. = Q x where is the A-dimensional sphere of radius p = \plE. 



10 



The billiard system has a natural equilibrium state given by the Liouville measure \i 
on A4, which is the product of a uniform measure on the polyhedron Q and a uniform 
(Lebesgue) measure on the sphere Sff , i.e. dfi = dqdp. The properties of billiard dynam- 
ics depend heavily on the curvature of the boundary dQ. In our case Q is a polyhedron, 
hence its boundary consists of flat sides with zero curvature. A prototype of such systems 
is billiard in a polygon. It is well known that (see, e.g., Q) 

Fact. For billiards in polygons and polyhedra (and hence, for our mechanical model of a 
piston in the ideal gas) all Lyapunov exponents vanish, and so does the Kolmogorov-Sinai 
entropy. 

Systems with zero Lyapunov exponents and zero entropy are not regarded as chaotic, 



but they still may be ergodic. In fact, billiards in generic polygons are ergodic ||KMS 
Moreover, for many nonergodic polygons, the phase space is foliated by invariant sub- 
surfaces on which the dynamics is ergodic. 

Even though there are no similar results, to our knowledge, for billiards in high- 
dimensional polyhedra, one can expect that they, too, have similar properties. That is, 
they are generically ergodic or become ergodic after trivial reductions. In our case, the 
billiard in Q is, perhaps, ergodic for typical values of M, or else the phase space is foliated 
by invariant submanifolds on which the dynamics is ergodic, and that those submanifolds 
fill M. pretty densely. In the latter case, one would hardly distinguish experimentally 
between such a nonergodic system and a truly ergodic one. 

Hence, we can assume that our system is ergodic or very close to ergodic in the 
above sense. Then almost every trajectory eventually behaves according to the invariant 
measure /i, independently of the initial state. In particular, for any initial gas density 
and velocity distribution (given by the function p(x, v), see Section ||) the hydrodynamic 
regime for a finite L is only valid on a finite interval of time - eventually the system will 
relax to equilibrium. We expect in fact that in terms of the "macroscopic" variables, say, 
the one particle distribution function, the system will relax to an effective equilibrium, 
as defined by ( |5.1| ) in terms of r eq , which is much smaller than the exponentially long 
time (in L) required for the ergodic theorem. So the real question is how does this time 
depend on L. According to our earlier discussion r c ~ logL and r eq ~ 

At equilibrium, the distribution of coordinates Xj and velocities Vi are determined 
by the Liouville measure /i, which is uniform in the phase space. Physically interesting 
(and only observable) are its marginal measures, i.e. projections, on lower-dimensional 
subspaces. The marginal measures of the velocities are normal (Gaussian) for large N 
(as iV — > oo). 

In particular, each individual velocity V{ converges in law to a Maxwellian (i.e., nor- 
mal) distribution with zero mean and variance 2E/N = const. The same holds for the 
"piston" component of the velocity, xq = Vy/M, hence the piston velocity V will be 
normally distributed with zero mean and standard deviation const/ yM =const/L, as 



L — > oo. In our case V has standard deviation J7/2A/L « 0.5/L. This conclusion agrees 
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well with our numerical datai. 

The equilibrium distribution of the piston coordinate X is also determined by the 
projection of the uniform measure dq on Q onto the xq axis. Before we do that, let us get 
rid of M in the definition of both Q and xq. A simple change of variable X = xo/y/M 
allows us to redefine Q by 

< X-n_, • • • , x-i < X < x\, . . . , xn + < L 

Furthermore, rescaling Y = X/L and yi = Xi/L gives a new, simpler, definition of Q: 

< y-N.i ■ ■ ■ , y-i < Y < 2/1, . . . , y N+ < 1 

"Integrating away" the variables yi yields the following equilibrium density for Y: 

f(Y) = cY N -(l-Y) N+ 

for < Y < 1, where c is the normalizing factor that can be computed explicitly. Put 
z = (Y — 0.5)^8/^, then the density of z is given asymptotically by f(z) ~ e~ z I 2 . 
Hence, Y is asymptotically gaussian with mean 0.5 and variance (4iV) _1 = (4L 3 ) -1 . 
Therefore, in equilibrium 



F-0.5 




and the probability of observing larger fluctuations is exponentially small in N. Note 
that fact is independent of the piston mass, as it has to be for an equilibrium (classical) 
system. This is consistent with our observation of the piston coordinate reaching a 
positive constant value, \Y — 0.5| ~ 0.1, during a short time interval < r < 20, since 
our initial conditions were selected according to the density function p(x, v) which has 
probability exponentially small in N . 

6 Remarks 

1. We checked the accuracy of our computer program in various ways. A simple one was 
based on the effect of round-offs on the energy of the system: the total energy was found 
to be practically constant over the whole period < r < 1000. A more sensitive test 
of the accuracy of the computations consists in using the time-reversal symmetry of the 
dynamics. Suppose at some time f > the velocities of all gas particles and the piston 
are reversed (vi — > —Vi). Then the system is supposed to trace back its past trajectory 
and arrive to its initial state with all velocities reversed at time 2f. We verified this 

5 It also allows us to estimate, in a peculiar way, the time of convergence to equilibrium, r oq , discussed 
earlier. Assume all the initial velocities are of order one, and note that the largest Maxwellian velocities 
arc of order y/N = L 3 / 2 . Since each collision with the piston adds 0(1/ L) to a particle's velocity, it 
takes ~ L 5 / 2 collisions to reach the maximum. An amazing agreement with the estimate r cq ~ L 2A7 
reported near Fig. 8. 
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property numerically for various f < 50 and always found that the system did repeat its 
past trajectory in the sense that the graphs of Y(r) and W{r) over the interval (f , 2r) 
looked like perfect mirror images of the corresponding graphs over the interval (0, r). 
Therefore, one can assume that round off errors remain negligibly small during times 
r < 50. However, this accuracy test failed for larger times, f > 100, indicating that the 
system then "loses memory" of its initial state due to round-offs. This does not seem to 
affect the overall picture, though. As yet another test, we carried out some computations 
with single rather than double precision, and found that the changes were little. Still, 
some estimates requiring long runs, such as the equilibrium time r eq on Fig. 8, might not 
be very accurate. 



2. We tested the dependence of the piston oscillations on the b factor involved in the 
formula Ml = bmL 2 . Remember that we set b = 2 in our main experiments. When we 
changed it to b = 20 (this made the piston 10 times heavier), then the oscillations started 
slightly later and their amplitude was slightly larger, but otherwise the picture was very 
much the same. When we changed b to 0.2 (and this made the piston 10 times lighter), 
then the oscillations started at about the same time as for 6 = 2, but they dampened 
somewhat faster. 

We also tried to change the piston mass Ml even more drastically. When we set it 
to L 3 (insteas of L 2 ), it appeared that oscillations started only after a very long initial 
quiescent period. But we did not examine this fact in detail, see ||GPL| . On the contrary, 



when we set the piston mass Ml to L (instead of L 2 ), large oscillations started very soon, 
but very quickly dampened and disappeared. 



3. We note that Eq. ( |3.3j ) describes the time evolution of a damped harmonic os- 
cillator. Accepting ( |3~3| ) over some time range, say, r G [30, 800] for L = 100, we can 
then look at the "inverse problem" of finding the effective spring constant and damping 
coefficient. 

Using the original variables, t and A^(t), we write 



Mr 



dt 2 



dX 

L +Kl(Xl-L/2)+til-^ 







which has the solution Xl — L/2 ~ e at with 



a 



2M T 



±i 







1 M l 


4K L M L _ 



.1/2 



This yields, to the lowest order in t^/AKlMl-, remembering that Ml = 2L 2 and that 
uj = 27r/r per , with r peT = t pcr /L 1.63, Kl ~ 8ti 2 /t 2 , i.e. the effective "restoring force" 
Kl is independent of L. On the other hand, the damping coefficient is tjl = 4L\l, were 
Xl was found experimentally to decrease as Xl ~ with 7 = 1.3, see Section 4. Hence 
t] L ~ L 1-7 = 0(L -0,3 ), i.e. the damping gets weaker as L — > 00. 
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7 Discussion 



We have presented here numerical results concerning the time evolution of a system with 
many degrees of freedom (up to 27 x 10 6 ) one of which, the position of the piston, plays 
a very special role. Starting with the particular initial particle distribution given by 
( [3.1D we found two striking features of the evolution: (1) the velocity distribution of the 
particles approaches a Maxwellian, i.e. the system goes toward thermal equilibrium and 
(2) the time evolution of the piston followed closely, after some initial period, that of a 
damped harmonic oscillator over an extended time interval, with initial oscillations as 
large as 1/10 of the system size. 

As already noted, we expect from general considerations |B| that (1) should be true 
for any initial density p(x,v) provided the number of particles N (~ L 3 ) is large enough 
and one waits long enough. But what about (2)? Clearly, if we choose for p(\v\) a 
Maxwellian, we expect only thermal fluctuations of order 0(L~ 3//2 ) in the position of the 
piston. This is indeed what we found numerically. For other, non-maxwellian, initial 
densities the question turns out to be far from trivial. We are currently working with 
more general initial densities and will report results in a separate paper ||CCLP|| . Below 
we outline our program and mention some preliminary findings. 

For simplicity, we assume that p(x, v) = p(\v\) and X(0) = L/2, V(0) = 0. In this 
case there is no apriori bias for the piston to move at all, and as already noted, the 
hydrodynamical (deterministic) equations, see ||CLS|| , predict that, in the limit L — > oo, 
the system would remain frozen in the initial state: Y(r) = Y(0) = 0.5 and W(r) = 
W(0) = for all t > (in the variables y and r, see Section 2). The density p(y,v,t) 
will also remain constant in time. But what about the particle system with a large but 
finite L? How will the piston and the gas behave while the particle velocities make their 
way to a Maxwellian? 

To answer this question we note that since the initial configuration of particles is 
generated randomly from a Poisson process with the density p(y, v,0) = p(\v\), the actual 
(empirical) density of the particles, such as the one shown on Fig. 6 at r = 0, does not 
exactly coincide with p(y, v, 0). Random fluctuations of the empirical density are typically 
of order 0(1/ L). Hence, the actual initial distribution of particles can be thought of as 
a small perturbation of the function p(\v\) and can be written as p(\v\) + epi(y, v) with 
e — l/L and some (random) function pi(y,v) of order one. 

Now, we conjecture that for large L the evolution of the particle system closely follows 
the solutions of the hydrodynamical equations derived in ||CLS| , |LPS|| with a perturbed 
initial density p(|i>|) +epi(y, v), rather than the stationary solution corresponding to the 
unperturbed density In particular, if the latter solution is unstable, then small 

perturbations grow exponentially in time (in the units of r = t/L, of course), hence 
the system can be destabilized in r ~ — loge = logL units of time. This would be in 
agreement with our analysis and numerical estimates in the end of Section 4. Hence, the 
instability of the hydrodynamical equations becomes a key issue. 

When we were finishing the present paper, we received a message from E. Caglioti 
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and E. Presutti who (a) proved that the hydrodynamical equations are stable when p(\v\) 
is monotonically nonincreasing in i.e. p'(M) — 0> an d (b) suggested that they might 
be unstable for our class on non-monotone p(|i>|). We checked the suggestion (b) for 
our particular density ( |3.1|) and found that it was indeed correct; we proved that small 
perturbations grow exponentially in r. Furthermore, when starting with a perturbed 
initial density with e = 0.01, we found numerically that the corresponding solution of 
the hydrodynamical equations resembled very well the evolution of the particle system 
described here, including large nearly harmonic oscillations of the piston during the 
interval 10 < r < 30. Further work in this direction is currently under way ||CCLP| . 



Conversely, when we simulated a particle dvnamics with a nonincreasing initial density 
p(H) the oscillations essentially disappeared!! On the other hand, the particle velocity 
distribution still approached a Maxwellian, albeit at a somewhat slower pace. 
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Figure 1: The piston coordinate Y as a function of time r. Here L = 100, iV_ = 500341, 
AT + = 499888. 
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Figure 2: The piston velocity W as a function of time r. The same run as in Fig. 1. 
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Figure 3: The piston coordinate Y during the intervals (30,35), (100,105), (250,255), 
and (900,905). The same run as in Fig. 1 and 2. 
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Figure 6: Six snapshots of the empirical gas density (in the x,v plane) at times r = 0, 
2.3, 4.2, 5.9, 7.4 and 18.6. 
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Figure 7: The value function of ln(W c L): experimental points and a linear fit. 
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Figure 8: The value lnr cq as a function of InL: experimental points and a linear fit. 




Figure 9: lnS* (thick line) and In S' (thin line) as functions of r. 
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